Conversation
|
Amazing. I can confirm it fixes not only the MWE I posted, but also the real code I extracted it from. Note that there is still a regression if the matrix is not wrapped in function barrier(point)
d2 = div(d, 2)
M = mat(point, T)
M[1:d2, d2+1:end] .= 0
M[d2+1:end, 1:d2] .= 0
return sum(real(eigvals(M)))
endIt doesn't matter to me, I never need a non-Hermitian matrix, I'm remarking just in case you wanted to know. |
| # analytic solutions | ||
| @test ForwardDiff.jacobian(eigvals_symreal, x) ≈ [0 0; 2 4] | ||
| @test ForwardDiff.jacobian(eigvals_hermreal, x) ≈ [0 0; 2 4] | ||
| if !(x isa StaticArrays.StaticArray) | ||
| @test ForwardDiff.jacobian(eigvals_symtridiag, x) ≈ [(1 - 3/sqrt(5))/2 (1 - 1/sqrt(5))/2 ; (1 + 3/sqrt(5))/2 (1 + 1/sqrt(5))/2] | ||
| end |
There was a problem hiding this comment.
I haven't looked very closely, but are there more tests somewhere that check the values?
In what's changed in this PR, I see one more line with Calculus.finite_difference_jacobian(ev, float.(x0)) where I think x0 = [1.0, 2.0].
There was a problem hiding this comment.
Yes, there are more in the lines below. Much more exhaustive than on master.
There was a problem hiding this comment.
Sorry that wasn't clear. I see the many tests of different paths against each other, which is good.
But what I meant is tests of any ForwardDiff path against something completely independent -- a known answer, or finite differences. To catch things like making the same logic error in writing both _eigvals_hermitian and _eigen_hermitian here.
There was a problem hiding this comment.
In the tests below, Jacobians of every test function with ForwardDiff are checked against finite differencing based Jacobians. The latter doesn't involve any ForwardDiff paths.
There was a problem hiding this comment.
Sorry to be dense but where are you looking? Besides line 293, as mentioned:
I see one more line with Calculus.finite_difference_jacobian(ev, float.(x0)) where I think x0 = [1.0, 2.0].
There was a problem hiding this comment.
It is line 293. Calculus doesn't support all vector types, hence the ForwardDiff results with x (the possibly static version of x0: lines 259 and 262) is compared with the finite differencing result using the regular array x0.
There was a problem hiding this comment.
Ok.
So the matrices being tested are these. Does it seem OK to test e.g. no negative eigenvalue for Symmetric{Real}, no zero eigenvalues for Hermitian{Complex}? Haven't thought hard but it just seemed a small sample.
julia> x
2-element Vector{Float64}:
1.0
2.0
julia> Symmetric(x*x') # also wrapped in Hermitian
2×2 Symmetric{Float64, Matrix{Float64}}:
1.0 2.0
2.0 4.0
julia> eigvals(ans)
2-element Vector{Float64}:
0.0
5.0
julia> Hermitian(complex.(x*x', x'*x))
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
1.0+0.0im 2.0+5.0im
2.0-5.0im 4.0+0.0im
julia> eigvals(ans)
2-element Vector{Float64}:
-3.0901699437494754
8.090169943749475There was a problem hiding this comment.
(That said, when I tried various tests locally, I did not manage to break this.)
|
The PR fixes #780 as well: #780 (comment) |
# Conflicts: # src/dual.jl
The matrices in the `eigen` testset are all of the form `x*x'` and hence symmetric, so `:U` and `:L` wrap the same matrix and the `uplo` handling these methods do is never exercised. Add a case whose raw storage is deliberately not symmetric. On master these fail: `eigvals(Symmetric(A, :L))` reads the values from the upper triangle -- `Symmetric(value.(parent(A)))` drops the `uplo` -- while reading the partials through the wrapper, so the two come from different matrices and the Jacobian is off by 0.198. The eigenvector case needs a matrix whose eigenvector direction actually moves: for `[x1 x2; x2 x1+x2]`, `tan(2θ) = 2x2/(a-d) = -2` is constant, so the derivative is identically zero and the test would pass either way. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`(Q' * ∂A) * Q` starts with `Adjoint * Symmetric`, which has no BLAS
specialization: it falls back to the generic `*(::AbstractMatrix,
::AbstractMatrix)`, allocating an extra n^2 temporary per partial and
skipping `symm`/`hemm` entirely (`mul!(C, Q', Symmetric(M))` is 1003 ns
against 169 ns for `mul!(C, Q', M)` at n = 10).
`∂A * Q` is `Symmetric * Matrix`, which does dispatch to `symm`, and the
remaining `Q' * _` is `Adjoint * Matrix`, i.e. `gemm`. Same result, one
fewer temporary, and the halved flop count now actually materializes.
This removes the one regression against master -- `eigvals(Symmetric)`
was 11.67 μs / 127 allocs / 58.5 KiB against master's 10.62 μs / 107 /
49.3 KiB, and is now 10.54 μs / 107 / 49.3 KiB -- and improves the rest
(gradient over a length-10 input):
master this branch
eigvals Symmetric 10.62 μs 49.3 KiB 10.54 μs 49.3 KiB
eigen Symmetric 25.54 μs 115.4 KiB 13.58 μs 68.6 KiB
eigvals Hermitian{C} 51.58 μs 544.8 KiB 25.75 μs 100.3 KiB
eigen SymTridiagonal 26.38 μs 89.4 KiB 14.00 μs 50.3 KiB
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Refreshed this branch: merged current Merged
|
master |
this branch | |
|---|---|---|
eigvals(Symmetric{<:Real}) |
10.62 μs, 107 allocs, 49.27 KiB | 10.54 μs, 107 allocs, 49.27 KiB |
eigen(Symmetric{<:Real}).values |
25.54 μs, 267 allocs, 115.41 KiB | 13.58 μs, 130 allocs, 68.56 KiB |
eigvals(Hermitian{<:Real}) |
MethodError |
10.83 μs, 107 allocs, 49.27 KiB |
eigen(Hermitian{<:Real}).values |
MethodError |
13.58 μs, 130 allocs, 68.56 KiB |
eigvals(Hermitian{<:Complex}) |
51.58 μs, 4133 allocs, 544.84 KiB | 25.75 μs, 133 allocs, 100.31 KiB |
eigen(Hermitian{<:Complex}).values |
MethodError |
33.38 μs, 156 allocs, 136.80 KiB |
eigvals(SymTridiagonal) |
11.23 μs, 163 allocs, 31.64 KiB | 10.56 μs, 143 allocs, 31.02 KiB |
eigen(SymTridiagonal).values |
26.38 μs, 380 allocs, 89.41 KiB | 14.00 μs, 166 allocs, 50.31 KiB |
No regressions left; the numbers in the PR description above predate current master and this change.
Two things deliberately not changed
_to_duals's complex branch computestuple.(parts...)twice, once for the real and once for the imaginary part. Hoisting it into a local looks like an obvious win but is a pessimization: inside the fused broadcast the tuples are never materialized, so hoisting forces an n²-element array of N-tuples into existence — 448 KiB against 240 KiB at n = 40, N = 8, to buy ~8% time. Left as is.- The repeated-eigenvalue case still divides by a vanishing eigenvalue gap and yields
NaN/Infeigenvector derivatives, exactly as onmaster(and as in ChainRules, which tracks it as an open TODO). That is Differentiate generaleigenandeigvals#788's territory, not this PR's.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #757 +/- ##
==========================================
+ Coverage 90.68% 91.23% +0.54%
==========================================
Files 11 11
Lines 1052 1072 +20
==========================================
+ Hits 954 978 +24
+ Misses 98 94 -4 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
3 participants
The PR adds missing definitions of
eigvalsandeigenforHermitian{<:Dual}and ofeigenforHermitian{<:Complex{<:Dual}}.Additionally, the PR adds more tests of
eigenandeigvals, unifies the existing implementations, and removes duplicate calculations in the definitions ofeigen.Fixes #756 without GenericLinearAlgebra.
For the example in #756, I get on ForwardDiff@0.10 without GenericLinearAlgebra:
And on ForwardDiff@0.10 with
import GenericLinearAlgebra: eigen:On this PR I get without GenericLinearAlgebra:
Some benchmarks against master:
master
This PR